## Use day of a week to define control group (Results are estimated by lfe package version 3.0-0)
## These results are reported in Table 7
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

#Using the last Saturday of December as the starting point - Date 69
data3$T[which(data3$Year==2017)]=data3$T[which(data3$Year==2017)]-1
data3$T[which(data3$Year==2016)]=data3$T[which(data3$Year==2016)]-2
data3$T[which(data3$Year==2015)]=data3$T[which(data3$Year==2015)]+3
data3$T[which(data3$Year==2014)]=data3$T[which(data3$Year==2014)]+2
data3$T[which(data3$Year==2013)]=data3$T[which(data3$Year==2013)]+1

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)
data_base1=data4
data_base2=data5
## Long sample
all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)
all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)
all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)
## Short sample
all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)
all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)
all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)
## Output results
stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C)
stargazer(all_SO2_short_C, all_NOX_short_C, all_AOD_short_C)
